{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fickian diffusion with variable diffusivity\n",
    "\n",
    "In this example, we will implement a Fickian diffusion algorithm on a `Cubic` network assumming a variable diffusivity, such that diffusivity is significantly higher at high concentrations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:29:58.891582Z",
     "iopub.status.busy": "2021-06-24T11:29:58.889931Z",
     "iopub.status.idle": "2021-06-24T11:29:59.668080Z",
     "shell.execute_reply": "2021-06-24T11:29:59.669192Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import openpnm as op\n",
    "%config InlineBackend.figure_formats = ['svg']\n",
    "import matplotlib.pyplot as plt\n",
    "np.random.seed(10)\n",
    "ws = op.Workspace()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating network, phase and geometry\n",
    "We create a `Cubic` network and add the geometry and phase objects. Note that the same procedure can be applied on extracted networks and 3D networks. Here for simplicity and visualization we chose a 2D cubic network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:29:59.680834Z",
     "iopub.status.busy": "2021-06-24T11:29:59.679408Z",
     "iopub.status.idle": "2021-06-24T11:29:59.815693Z",
     "shell.execute_reply": "2021-06-24T11:29:59.814293Z"
    }
   },
   "outputs": [],
   "source": [
    "net = op.network.Cubic([1, 15, 15], spacing=1e-5)\n",
    "geom = op.geometry.SpheresAndCylinders(network=net, pores=net.Ps, throats=net.Ts)\n",
    "air = op.phases.Water(network=net)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding the diffusivity model\n",
    "We assume that the diffusivity of the phase air is not constant and is a function of the concentration of the species. Therefore, we define a model to calculate the `pore.diffusivity` at each concentration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:29:59.826173Z",
     "iopub.status.busy": "2021-06-24T11:29:59.824760Z",
     "iopub.status.idle": "2021-06-24T11:29:59.828764Z",
     "shell.execute_reply": "2021-06-24T11:29:59.829915Z"
    }
   },
   "outputs": [],
   "source": [
    "def variable_diffusivity(target, c_ref, pore_concentration=\"pore.concentration\"):\n",
    "    X = target[pore_concentration]\n",
    "    val = 1e-9 * (1 + 5 * (X / c_ref) ** 2)\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:29:59.845208Z",
     "iopub.status.busy": "2021-06-24T11:29:59.843774Z",
     "iopub.status.idle": "2021-06-24T11:29:59.881701Z",
     "shell.execute_reply": "2021-06-24T11:29:59.882863Z"
    }
   },
   "outputs": [],
   "source": [
    "air.add_model(propname=\"pore.diffusivity\", model=variable_diffusivity, c_ref=10.0,\n",
    "              regen_mode=\"deferred\")\n",
    "phys = op.physics.Standard(network=net, geometry=geom, phase=air)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Performing Fickian diffusion\n",
    "We then add the the FickianDiffusion algorithm to our simulation and define \"Dirichlet\" boundary conditions assigned to the  pores that are on the top and the bottom face of the network. Note that values=10.0 for the first boundary condition means that the pore concentration of each pore (that is on the `top` boundary face) is equal to 10.0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:29:59.890707Z",
     "iopub.status.busy": "2021-06-24T11:29:59.889239Z",
     "iopub.status.idle": "2021-06-24T11:30:00.254396Z",
     "shell.execute_reply": "2021-06-24T11:30:00.253808Z"
    }
   },
   "outputs": [],
   "source": [
    "fd = op.algorithms.FickianDiffusion(network=net, phase=air)\n",
    "fd.set_value_BC(pores=net.pores('top'), values=10.0)\n",
    "fd.set_value_BC(pores=net.pores('bottom'), values=0.0)\n",
    "fd.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing\n",
    "\n",
    "We can visualize the average concentration distribution along the x direction (from the left boundary face towards the right boundary face)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:00.277593Z",
     "iopub.status.busy": "2021-06-24T11:30:00.273046Z",
     "iopub.status.idle": "2021-06-24T11:30:00.394701Z",
     "shell.execute_reply": "2021-06-24T11:30:00.395195Z"
    }
   },
   "outputs": [],
   "source": [
    "air.update(fd.results())\n",
    "c=air['pore.concentration']\n",
    "fig, ax = plt.subplots()\n",
    "c_xz_avg = c.reshape(op.topotools.get_shape(net)).mean(axis=(0, 1))\n",
    "ax.plot(c_xz_avg, \"ko-\")\n",
    "ax.set_xlabel(r\"z\")\n",
    "ax.set_ylabel(r\"c\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The concentration of species throughout the network is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:00.418268Z",
     "iopub.status.busy": "2021-06-24T11:30:00.411358Z",
     "iopub.status.idle": "2021-06-24T11:30:00.580330Z",
     "shell.execute_reply": "2021-06-24T11:30:00.579758Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x2b9ce6b4a60>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c2d = c.reshape(op.topotools.get_shape(net)).squeeze()\n",
    "plt.imshow(np.rot90(c2d))\n",
    "plt.colorbar()"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3 (Spyder)",
   "language": "python3",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
